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Abstract 

We present a field formulation for defects that draws from the classical representation of 
the cores as force dipoles. We write these dipoles as singular distributions. Exploiting the key 
insight that the variational setting is the only appropriate one for the theory of distributions, 
we arrive at universally applicable weak forms for defects in nonlinear elasticity. Remarkably, 
the standard, Galerkin finite element method yields numerical solutions for the elastic fields of 
defects that, when parameterized suitably, match very well with classical, linearized elasticity 
solutions. The true potential of our approach, however, lies in its easy extension to generate 
solutions to elastic fields of defects in the regime of nonlinear elasticity, and even more notably 
for Toupin’s theory of gradient elasticity at finite strains (Arch. Rat. Mech. Anal., 11 , 385, 
1962). In computing these solutions we adopt recent numerical work on an isogeometric analytic 
framework that enabled the first three-dimensional solutions to general boundary value problems 
of Toupin’s theory (Rudraraju et al. Comp. Meth. App. Mech. Engr ., 278 , 705, 2014). We 
first present exhaustive solutions to point defects, edge and screw dislocations, and a study on 
the energetics of interacting dislocations. Then, to demonstrate the generality and potential 
of our treatment, we apply it to other complex dislocation configurations, including loops and 
low-angle grain boundaries. 

1 Introduction 

The elastic fields of crystal defects have proven difficult to represent. As is well appreciated, 
the large distortion in the core places the elasticity problem in the nonlinear (i.e., finite strain) 
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regime. Furthermore, due to the very non-uniform distortion in the core, measures that go be¬ 


yond the deformation gradient exert a magnified influence. Generalized (Cosserat and Cosserat 


1909), higher-order (Toupin, 1962) and nonlocal (Eringen and Edelen, 1972) theories of elas¬ 


ticity that also preserve geometric nonlinearities therefore provide the appropriate treatment. 
Here, we work with Toupin’s theory of gradient elasticity at finite strains (Toupin 1962) moti¬ 
vated primarily by the emergence of strain gradients as second-order terms in the Taylor series 
expansion of the elastic free energy density (Rudraraju et ah, 2014). 

A number of different approaches have been adopted to obtain analytic solutions to clas¬ 
sical (i.e., without higher-order measures of deformation), nonlinearly elastic, boundary value 
problems of defect configurations. Rosakis and Rosakis (1988) developed anti-plane shear field 
solutions to screw dislocations in nonlinear, homogeneous, isotropic, incompressible solids. Their 
work defines a stress-based measure of nonlinearity, motivated by the differences in total stored 
energy that arise, depending on the chosen strain energy function. ? developed a treatment of 
distributed dislocation and disclination densities at finite strains by invoking internal degrees of 
freedom and couple stresses. Analytic solutions were obtained in the setting of plane elasticity. 
Working with distributed defects within the setting of differential geometry, |Yavari and Goriely| 
(2012) described analytic solutions to a number of screw and edge dislocation configurations, 
finding, in some cases, that the singularity depends on the material model. Later, these au¬ 
thors extended their analysis to point defects (?). By representing the cores with a continuous 
distribution of defects, the authors avoided the singularities that emerge in the classical lin¬ 
ear solutions. Using a continuum theory of geometrically necessary dislocations, ? obtained 
analytic solutions to screw dislocations in a neo-Hookean solid. Working with the continuum 
theory of dislocations based on a weak definition of the dislocation density tensor, |Acharya] 
(2004) subsequently developed singularity-free field solutions in a partial differential equation 
framework. Subsequently, the author presented a statistical mechanics-based thermodynamic 
theory of dislocations, also using the dislocation density tensor (Acharya, 2011). A number of 
other significant works have been cited in the above references; we are constrained only by space 
from listing more of them here. 

Significant as the above works are, there remains room for a numerical treatment of the 
problem of defect fields in classical, nonlinear elasticity, especially in finite bodies, and for more 
complex defect configurations. To develop such a framework is a goal of our work, but more 
crucially, we seek a framework that can be extended to gradient elasticity at finite strains as a 
means to regularize the core fields. This is not a straightforward task as we aim to argue. 

Such solutions that do exist of higher-order elasticity theories applied to defect structures 
are restricted to a linearized treatment. Typical, are those based on Mindlin’s formulation of 
linearized gradient elasticity (Mindlin 1964 ), further reduced to lower dimensions and simplified 
boundary value problems. There are a number of such treatments, of which we point the reader 
to Gutkin and Aifantis (1999); Lazar and Maugin (2005). Also see Kessel (1970); Lazar et al. 


(2005, 2006) for applications of generalized theories of elasticity to defect fields in the linear 
regime. Higher-order elasticity theories, because of their complexity, have continued to resist 
general solution in three dimensions. This challenge has kept even numerical, field solutions out 
of reach to higher-order, finite strain elasticity treatments of defects in general, three-dimensional 
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boundary value problems. 

The present work relies on a numerical framework, developed in Rudraraju et al. (2014), 
which enables solutions of Toupin’s gradient elasticity theory at finite strain. It is based on 
isogeometric analysis (Hughes et ah, 2005 Cottrell et al. 2009) and, in particular, exploits the 
ease of developing spline basis functions with arbitrary degree of continuity in this framework. 
This is critical because, in the weak form of strain gradient elasticity, second order spatial deriva¬ 
tives appear on the trial solutions as well as the weighting functions, requiring functions that 
he in ^ 2 . This requirement is satisfied by C 1 -continuous spline functions. The framework in 


Rudraraju et al. (2014) demonstrated the first three-dimensional solutions to general boundary 


value problems of Toupin’s theory of gradient elasticity at finite strains. We seek here to exploit 
this numerical treatment and extend the catalogue of three-dimensional solutions of Toupin’s 
theory to include defect fields. 

Classical descriptions of defect fields, such as the Volterra dislocation (Volterra 1 1907), are 
based on a representation of the defect field away from the core (the far-held) that incorporates 
a displacement discontinuity at the core. Absent a direct model of the core, the elastic fields are 
not accurate near the core. A different approach is to represent the core by dipole arrangements 
of forces, which model the tractions experienced by the atoms immediately surrounding the 
core. Such a model has been employed for point defects to derive the Green’s function solution 
(Hirth and Lothe, 1982). For dislocations, the potential-based methods of Eshelby et al. (1953) 


have been applied to determine the linear elastic fields in Gehlen et al. (1972) and Hirth and 


Lothe (1973) by modelling the force dipole as an ellipsoidal center of expansion. In Sinclair 


et al. (1978) lattice Green’s functions were used in conjunction with the analytic expressions 


of Eshelby and co-workers to couple the far-held elasticity solution with the core held obtained 
from molecular statics. More recently, this class of approaches was extended to compute the 
elastic energies of dislocations by parameterizing the linear elastic helds against ab initio cal¬ 
culations (Clouet 2011 Clouet et al., 2011). These works rely on analytic expressions or series 
expansions of the elastic helds. Here, we aim to incorporate dipole force representations within 
large scale, mesh-based numerical methods such as hnite element or isogeometric methods. The 
key insight that allows us to realize this goal is that dipole force arrangements have a singular 
distributional character. In fact, dipole force arrangements give rise to singular dipole distri¬ 
butions. Furthermore, the variational setting is the only correct one for singular distributions. 
Therefore, variationally based numerical techniques are particularly well-suited to compute the 
helds resulting from the singular dipole distributions. This holds for hnite element and isogeo¬ 
metric methods, with the only distinguishing feature being that isogeometric methods ease the 
representation of C 1 -functions for strain gradient elasticity. 

We hrst review Toupin’s theory of gradient elasticity at hnite strain in Section [2j This is 
followed by a derivation of the singular distributional representation for dipole arrangements of 
forces in Section [3] The numerical framework, consisting of isogeometric analysis and including 
quadrature rules, appears in Section [4] An extensive set of numerical results for the elastic 
helds of point and line defects is presented in Section [5] We conclude by placing our work in 
perspective in Section [6] 
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2 Toupin’s theory of strain gradient elasticity at finite 
strains 


2.1 Weak and strong forms 

Our treatment is posed in the Cartesian coordinate system, with basis vectors e$, i = 1,... 3, 
ei • Bj = Sij . The reference configuration, its boundary and the surface normal at any boundary 
point are denoted by fio, SUo and IV, respectively, with |JV| = 1. The corresponding entities 
in the current configuration are denoted by U, dd and n, respectively. We work mostly with 
coordinate notation. Upper case subscript indices are used to denote the components of vectors 
and tensors in the reference configuration and lower case subscript indices are reserved for those 
in the current configuration. Working in the reference configuration, we consider the boundary 
to be the union of a finite number of smooth surfaces Tq, smooth edges Tq and corners Sq: 
dflo = To U To U So, for full generality. For functions defined on 9Uo, when necessary, the 
gradient operator is decomposed into the normal gradient operator D and the surface gradient 
operator Dx, 


= DtJjNk + D K !p 


where D^Nk = V’ jNjNk and Dk^> = 


(1) 


A material point is denoted by X E do. The deformation map between do and d is given 
by ip(X,t) = X + u = cc, where u is the displacement field. The deformation gradient is 
F = dtp/dX = 1 + du/dX , which in coordinate notation is expressed as Fu = dpi/dXj = 
Sij + dui/dXj. The Green-Lagrange strain tensor is given in coordinate notation by Eij = 
\{F kI F kJ - Sij). 

The strain energy density function is W(E, Gradi£). We recall that the dependence on E 
and Grad E renders W a frame invariant elastic free energy density function for materials of 


grade two (Toupin 1962). Constitutive relations follow for the first Piola-Kirchhoff stress and 
the higher-order stress tensors: 


Pu = 


dW 

dF~j 


( 2 ) 


BiJK — 


dW 

mJ, 


K 


(3) 


We note that instead of Equations © and © the single stress tensor that combines them 


could be used as in Toupin (1964). Our treatment (Rudraraju et ah, 2014, 2015) has been 


based on Toupin (1962), and in a future communication we will consider a comparison of the 


two approaches for the representation of defects. One obtains the same governing equations, 


constitutive prescriptions and jump conditions across internal interfaces as in Toupin (1964) 
with a strain energy density function of the form W(E, Gradl£), but no higher-order stress in 
the theory, and a global statement of the second law. See |Acharya and Fressegeas (2015) in this 
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regard, where the authors recover the formulation of Toupin (1964) with couple stress, but not 
higher-order stress. 

In the most general case, we have a body force distribution /, a surface traction T, a 
surface moment M and a line force L. For i = 1,2,3 denoting the Cartesian coordinates, 
the smooth surfaces of the boundary are decomposed as To = Tq* U = T£™ U , and the 
smooth edges of the boundary are decomposed as T 0 = YC U Y^. Here, Dirichlet boundary 
subsets are identified by superscripts u , m and l and Neumann boundary subsets are identified 
by superscripts T, M and L. 

We begin with the weak form of the problem. We seek a displacement field of the form 


Ui G with Ui = Ui on Tq*; U{ = on Yq* ; 


Dui = rni ; on TJJJ 


(4) 


such that for all variations of the form 


G 'V , with Wi = 0; on T 


U Yq*, Dwi = 0; on 


(5) 


the following equation holds: 

f (wi,jPu + w it j K BijK ) dV— [ Wifi dV- [ WiTi dS- [ DwiMi dS- [ WiL t d C = 0. 
Jn 0 Jn 0 Jrl J r" J 

(6) 

The problem has a fourth-order character, which resides in products of Buk and each 

of which involves second-order spatial gradients. Standard variational arguments lead to the 
strong form of the problem: 


PiJ,J - B iJK ,JK + fi 

= 0 

in f2o 

Ui 

= Ui 

on 

PijNj - DB uk N k Nj - 2 Dj(B iJK )N K - B iJK DjN K + ( b L L NjN K - b JK )B iJK 

= Ti 

on 

Dui 

— fhi 

on 

BukNjNk 

= Mi 

on r£f 

Ui 

= U 

on 

IN V jN k B iJK \ 

= Li 

on Tgi 


(7) 

Here, bu = —DiNj = —DjNi are components of the second fundamental form of the 
smooth parts of the boundary and N r = S x IV, where S is the unit tangent to the curve 
Yo (Toupin 1962). If Yo is a curve separating smooth surfaces Tq" C To and T^" C To, with 
N r+ being the unit outward normal to Yq from Tq" and N r being the unit outward normal 
to Y 0 from Tq we define \N^NkBuk\ := Nj + NkB^jk + Nj NkBuk- The (nonlinear) 
fourth-order nature of the governing partial differential equation above is now visible in the 
term Buk,jk , which introduces F a B,cjK via Equation Q. The Dirichlet boundary condition 
in 02 has the same form as for conventional elasticity. However, its dual Neumann boundary 
condition 03 is notably more complex than its conventional counterpart, which would have 
only the first term on the left hand-side. Equation 04 is the higher-order Dirichlet boundary 
condition applied to the normal gradient of the displacement field, and Equation Qs is the 
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higher-order Neumann boundary condition on the higher-order stress, B. Adopting the physical 


interpretation of B as a couple stress (Toupin 1962), the homogeneous form of this boundary 


condition, if extended to the atomic scale, states that there is no boundary mechanism to impose 
a generalized moment across atomic bonds. Finally, Equation 06 is the Dirichlet boundary 
condition on the smooth edges of the boundary and Equation 07 is its conjugate Neumann 


boundary condition. Following Toupin (1962), the homogeneous form of this condition requires 


that there be no discontinuity in the higher order (couple) stress traction across a smooth edge 
Tq in the absence of a balancing line traction along Tq. In Rudraraju et al. (2014) we have 
detailed the variational treatment leading to Equations 0 and 0, as well as the corresponding 
statements in the current configuration. 


3 Representation of defects as force dipole distributions 

3.1 Point defects 

In fl 0 , consider a discrete dipole formed of force vectors — R and R located at distinct points X[ 
and X' 2 , respectively. The force field is a singular distribution. While the theory of distributions, 
rigorously applied, holds in a variational setting, it permits a formal representation as classical 
functions (Stakgold 1979| ). This allows a manipulation of the force fields as follows: 


/(X; X[, X ' 2 ) = RS 3 (X; X') - RS 3 (X; X[), 


(8) 


where S 3 (X;X / ) is the three-dimensional Dirac delta distribution. Expanding it in a Taylor’s 
series around Y' = (X[ + X r 2 )/2. 


RS 3 (X ; X\) = R[5 3 (X ; Y')\ + R ( ^7<5 3 (X; X') 


RS 3 (X; X') = R[6 3 (X; Y')] + R ( ^b<J 3 (X; X') 


Y 


(xi-v'j + oflxi-y'i 2 ) (9) 


Y 


• (X' -Y') + o(\x' 2 -V| 2 ), ( 10 ) 


To first order, therefore, the force field is a dipole distribution 

d 


f(X;Y')=R( M 7 S 3 (X;X) 


Y 


■ (*2 - *!)• 


Letting £ = X' 2 — X\ for brevity, we have, in coordinate notation, 

dS 3 (X;X') 


fi = Ri 


dX' 


Y 




We define the dipole tensor, 

D = R®£, 

and use cW 3 (X; X')/dX = — <9<5 3 (X; X')/dX' to write the force distribution, 

cM 3 (X; Y') 


f(X:Y') = -D- 


dX 


( 11 ) 


( 12 ) 


(13) 


(14) 


6 






























Figure 1: A dipole of forces representing a substitional point defect. 
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A direct application of the above development is the representation of a point defect as a center 
of expansion or contraction. In this case, the dipole is a diagonal tensor (Figure [l]). Using the 
ortho normal Euclidean basis {ei,e 2 ,es}, we have 

3 

D = ^2 R i£i e i® e u ( 15 ) 

i,I=1 

where Ri ^ R2 ^ R3 and (1 ^ (2 7 ^ £3 to model anisotropic point defects. 

3.2 Line defects 



Figure 2: 


An edge dislocation represented as an interfacial distribution of force dipoles. 


We retrace the above development from Equation (| 8 | through (14), but with interfacial 
force densities —R 1 and R 1 that are uniform over parallel planes r[ and respectively. In 




















this case, we arrive at relations with the same form as Equations (14) and (15). However, the 


one-dimensional Dirac-delta distribution, S 1 (X' 1 X') replaces the three-dimensional Dirac-delta 
distribution in the expression for the force distribution, 


f{X;Y') = -D 1 


dX 


(16) 


The resulting interfacial distribution of force dipoles is shown in Figure [2] for an edge dislocation, 
where the dipole tensor is defined as 

3 

£>*= £) (17) 


i,I=1 


Here, £ denotes the vector pointing from T[ toward T' 2 , R 1 = for edge dislocations, and 
R 1 • £ = 0 for screw dislocations. 


3.3 Linearized elasticity solutions for the dipole tensor 

In the limit of infinitesimal elasticity it can be shown that the dipole tensor of a point defect in 


an infinite domain is (Garikipati et al., 2006) 


Dij = CijklVk L , 


(18) 


where V T is the relaxation volume tensor of the point defect. Note that coordinate notation in 
this section uses uppercase superscript indices because of the coincidence of and Qq in the 
infinitesimal limit. 

A related result is possible for line defects by invoking Volterra’s linearized elasticity solution 
for the displacement field of a dislocation in an infinite domain: 

dG M K {X:X r ) t 


■ m =-L 


-CuKL^ibjdS 


(19) 


lM ^ 

where Gmk = Gkm is the infinite space Green’s function for elasticity, and gives the displace¬ 
ment in the cm direction for a unit force in the direction. The surface of integration is 
the half plane of the line defect, N is the unit normal to T' and b is the Burgers vector of the 
dislocation. The definition of the Green’s function implies that, also can be written using 
the force distribution introduced in ©. 


h M 


(X') = - f G MK (X;X')f K (X,X")dV 
Jn 


( 20 ) 


Using Equations (16) and and a standard result on the gradient of the Dirac-delta distri¬ 
bution 

«m( X') = - [ Gmk(X; X'){-D' kl ^% )dV 
Jn 


f dGMKj^JG) 

Jn 0X L 

[ 9G M k(X; X 1 ) 

Jr" 


dX L 


{D^ l 6\X-,X")) dU 


dX L 


X=X" 


D [ KL dS r ‘ 


( 21 ) 
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Comparing Equations (19) and (21) we obtain 


D'kl = CiJKL^jbj 


( 22 ) 


Using Cijkl = ^ij^kl + n(ftiK&jL + 5il3jk) in the isotropic case, and b = (3N for an edge 
dislocation, we have 


D% l = MklNA + v(N K b L + 7 V l 6 x ). 


(23) 


Using b ■ N = 0 for a screw dislocation, leads to 

Drl = n(N K b L + N L b K ). 


(24) 


Equations (18), (23) and (24) allow a direct comparison of our approach with linearized 


elastic fields for the corresponding defects. In the remainder of this communication, we will 
continue to use these relations as approximations to the respective point and line defect dipole 
tensors outside the linearized elastic regime of infinitesimal strain; i.e., for nonlinear elasticity 
with and without gradient effects at finite strains. In drawing conclusions at the end of the 
manuscript, we outline our approaches for better estimates of these defect dipole tensors. 

3.4 A field formulation for defects in gradient elasticity at finite strains 

On substituting Equation (14) and (16) in the weak form ©, we have, for point defects: 

dS 3 (X; Y') 


r r 

/ + BiJKWi,JK ) dV+ Wi {X)Du- 

J Q n J On 


dx. 


dF- 


f WiTi dS — [ DwiMi dS — f WiLi d C = 0, 

J r T Jv M J r L 


(25) 


and for line defects, 


J dXj dF 


/ (Puwi,j + Bukw^jk) dF+ f w i (X)D i d6l ( X ’ Y ^ 

J O 0 J Oo 

[ d S - [ Dw^i d S - [ w^i d C = 0 . 

J r T J r M J t 


(26) 


On again using the result for the gradient of distributions, and the fundamental definition of 
the Dirac-delta distribution these simplify to 

/ + Bukw^jk) dU - w i: j(Y')Dij - 

J o 0 

[ d S - [ D Wi Mi dS - / w^i d C = 0, (27) 

J r T . Jt m . J t l . 
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for point defects, and 


/ {PijWi,J + Bij K Wi,jK ) dV — / Wi t jD\j 

Jn o Jr' 


dS- 


[ WiTi dS - [ 
Jr T . Jr* 

O' 1 


DwjMj d S — 


l 

Jr 1 


WjLj dC = 0 


(28) 


for line defects, with Du given by (18), and D\j given by(23) or (24) 


It bears emphasizing that the rather transparent form of (27) and (28) is a direct consequence 


of the variational nature of the theory of distributions. The singular, Dirac-delta and dipole 
distributions are meaningful only in the variational setting. With a weak form at hand, their 
actions are simply transferred to the variations—or the test functions, w. We proceed to show 
that this leads to a surprisingly effective representation of point and line defect fields using 
variationally based numerical methods. 

4 Numerical treatment 


4.1 Galerkin formulation 

As always, the Galerkin weak form is obtained by restriction to finite dimensional functions (•) h : 

Find G y h C 5?, where 5? h = {u^ G JP 2 (Qq) \ = tq on T^, Du% = fhi on u\ = Z^on T^}, 

such that V G Y h C Y, where Y h = {w 1 ^ G Jrf? 2 (Q o) \ w i = 0 on T^UY^, Dw f = 0 on T™} 


/ W 

J On 


J W i,J 




iJK w i,JK ) 


dV- 


W i,jDiJ 


Y' 


(29) 


(30) 


- f w’lTi d S- [ Dw^Mi AS - f w'lU d C = 0, 
Jr T . Jr M J r L . 

QX QZ QZ 

for point defect fields, and 

f (Pfjwt + B? JK ™t K ) d V- [ w^jD\j as 

Jrtn Jr' 


- f as - [ Dv/'M, as- [ v4 1 l, ac = o, 

Jr T . Jr M Jr L . 

QZ QZ QZ 

for line defect fields. 

Our adoption of Toupin’s theory of gradient elasticity at finite strains is motivated by an 
interest in obtaining fields that are accurate at large strains, while remaining singularity-free 
through the core. As is well-appreciated, the second-order gradients in the weak form re¬ 
quire the solutions to lie in a more restrictive condition than the formulation of finite 

strain elasticity for materials of grade one, where the solutions are drawn from the larger space 
D The variations, w h and trial solutions u h are defined component-wise 

using a finite number of basis functions, 

n b n b 


W 


= ca N a , u h * ^2 daNa 


(31) 


a= 1 


a=l 


li 














where n b is the dimensionality of the function spaces and V h , and N a represents the basis 
functions. Since 5? h C ^ 2 , C° basis functions do not provide the degree of regularity demanded 
by the problem; however, it suffices to consider C 1 basis functions in . One possibility is the 
use of C 1 Hermite elements as in Papanicolopulos et al. (2009). Alternately, one could invoke 
the class of continuous/discontinuous Galerkin methods Engel et al. (2002); Wells et al. (2004); 
Molari et al.| ( | 20 06| ; | Wells et al.| ( |2006[ ), in which the displacement field is C°-continuous, but 


the strains are discontinuous across element interfaces. A mixed formulation of finite strain 
gradient elasticity could be constructed by introducing an independent kinematic field for the 
deformation gradient or another strain measure. These last two approaches, however, incur 
additional stability requirements. We prefer to avoid the complexities of Hermite elements in 
three dimensions, and seek to circumvent the challenges posed by discontinuous Galerkin meth¬ 
ods and mixed formulations by turning to Isogeometric Analysis introduced by |Hughes et al. 
(2005). Also see Cottrell et al. (2009) for details. 


4.1.1 Isogeometric Analysis 


As is now well-appreciated in the computational mechanics community, Isogeomeric Analysis 
(IGA) is a mesh-based numerical method with NURBS (Non-Uniform Rational B-Splines) basis 
functions. The NURBS basis leads to many desirable properties, chief among them being the 
exact representation of the problem geometry. Like the Lagrange polynomial basis functions 
traditionally used in the Finite Element Method (FEM), the NURBS basis functions are par¬ 
titions of unity with compact support, satisfy affine covariance (i.e an affine transformation of 
the basis is obtained by the affine transformation of its nodes/control points) and support an 
isoparametric formulation, thereby making them suitable for a Galerkin framework. They enjoy 
advantages over Lagrange polynomial basis functions in being able to ensure C n -continuity, in 
possessing the positive basis and convex hull properties, and being variation diminishing. A de¬ 
tailed discussion of the NURBS basis and IGA is beyond the scope of this article and interested 


readers are referred to Cottrell et al. (2009|). However, we briefly present the construction of 
the basis functions. 


The building blocks of the NURBS basis functions are univariate B-spline functions that 
are defined as follows: Consider two positive integers p and n, and a non-decreasing sequence 
of values x = [£i? £ 2 ? • •••, Cn+p+i] 5 where p is the polynomial order, n is the number of basis 
functions, the ^ are coordinates in the parametric space referred to as knots (equivalent to 
nodes in FEM) and x the knot vector. The B-spline basis functions Bi jP (£) are defined 
starting with the zeroth order basis functions 


* 5(0 = 


1 if &<£<&+i, 
0 otherwise 


and using the Cox-de Boor recursive formula for p > 1 (Piegl and Tiller} |199 7|) 


m) 


Ci+p C 


+ A +P +1 / B;t 1(0 


(32) 


(33) 
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The knot vector divides the parametric space into intervals referred to as knot spans (equiv¬ 
alent to elements in FEM). A B-spline basis function is C °°-continuous inside knot spans and 
C p ~ l -continuous at the knots. If an interior knot value repeats, it is referred to as a multiple 
knot. At a knot of multiplicity &, the continuity is C p ~ k . Now, using a quadratic B-spline basis 
(Figure ([3|), a C 1 -continuous one dimensional NURBS basis can be constructed. 0 



Figure 3: A quadratic B-spline basis constructed from the knot vector x 

[0,0, 0,1/6,1/3,1/2, 2/3,5/6,1,1,1]. 






(34) 


where Wi are the weights associated with each of the B-spline functions. In higher dimen¬ 
sions, NURBS basis functions are constructed as a tensor product of the one dimensional basis 
functions: 




£”=i E”=i B p(0 B i(v)wij 


(2D) (35) 




B^Bj^B^ Qwijk 

£"=i £"= 2 1 £ft=i B^)Bi( V )Bl(()w ijk 


(3D) (36) 


4.2 Numerical integration of singular force distributions 


From the theory of distributions, the forcing term is applied at the point defect in equation 
(29) and along the plane of the dislocation in equation (30). Special quadrature points must be 
introduced to numerically integrate these terms. For the point defect in Equation (29) this is 
accomplished with a single quadrature point as shown in Figure @ and Equation ( |37| ). 


w 




Y 


n b 

a= 1 


Y' 


.Dt 


(37) 


1 The boundary value problems that follow in Section [H] consider only simple geometries. For this reason, we 
have used the simpler B-spline basis functions instead of the NURBS basis. However, we have included the latter 
in this discussion for the sake of completeness, noting that the numerical formulation as presented is valid for any 
single-patch NURBS geometry. 
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Figure 4: Quadrature for a point defect. 
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The scheme for integration along the dislocation plane in Equation (30) is similar to that for 
Neumann boundary conditions. An internal surface of quadrature points is introduced on the 
plane as shown in Figure [5] and Equation (38). Two-dimensional Gaussian quadrature is found 
to be sufficient since the B-spline-derived functions are polynomials. The forcing term is, 



Figure 5: Quadrature points on a dislocation plane. 


„ n q n h n h n q 

/ KjDijds = e E = E < E 
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where n q is the number of quadrature points on the plane T', and £ q and w q are coordinates 
and weights for the corresponding quadrature points. 

5 Numerical results 

For the full generality of Toupin’s theory, we consider an elastic free energy density function W, 
that incorporates gradient effects at finite strain. As is well-known, material frame invariance is 
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guaranteed by requiring W to be a function of the Green-Lagrange strain tensor, E = |( F T F — 
1), and its gradient. We choose a simple extension of the St. Venant-Kirchhoff function by a 
quadratic term in Gradi£, and write it in coordinate notation: 

W = (E, Grad E) = — ( Eaa ) 2 + H^EabEab) + 1 2 Eab,cEab,c , (39) 

where l is a gradient length scale parameter. The Lame parameters have been normalized so 
that A = 1 and fi = 1. All computations were carried out on the unit cube (0, l) 3 , relative to 
which the magnitude of the Burgers vector 5,, and gradient length scale l have been normalized. 

To validate the dipole representation of defects, we first simplify our model to linearized 
elasticity (infinitesimal strain) and ignore the gradient effects. This enables a comparison be¬ 
tween the classical, analytical solutions in this regime and numerical solutions for the point 
defect, edge dislocation and screw dislocation. We also include in this study, comparisons with 
solutions obtained with gradient elasticity at finite strains, which demonstrate the expected reg¬ 
ularization of displacement and stress fields. Following this establishment of the fundamental 
solution characteristics, we present a study of the strain energy of a single dislocation, and of 
two interacting dislocations in the setting of classical, linearized elasticity as well as gradient 
elasticity at finite strains. Finally, to demonstrate the potential for extension to more complex 
defect configurations, we compute the solutions for an edge dislocation near a traction-free sur¬ 
face, a non-planar dislocation loop, and a grain boundary represented by an arrangement of 
edge dislocations. 


5.1 Comparison of analytic and numerical solutions for point defects 

In the regime of linearized elasticity, the analytical solution to the displacement field around a 


point defect in an infinite medium is (Hirth and Lothe 1982) 


.an 
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(40) 


where the dipole tensor is D — D 1 , with the normalization D — 1 x 10 . To relate D to the 


relaxation volume tensor, we have from (18) in the isotropic case, 


D IJ = \6 I jVk K + AVIj + VJ I ) 


(41) 


In this 


Vh = 


1 


3A 2 /j, 


DSu 


(42) 


All our numerical solutions are over a unit cube = (0, l) 3 . For the point defect, we consider 
an interstitial located at X' = {0.5, 0.5, 0.5} T . For consistency, we apply the analytical solution 


on the boundary surfaces Tq as Dirichlet boundary conditions. The weak form (27) is then 
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re-written as: 


I wi'jvudV - w iti 


(43) 


X 


,D = 0 


u = u an on r 0 . (44) 

Numerical solutions are obtained using isogeometric analysis as described in Section]!} Figure 
[6] is the deformed configuration of the cube around the interstitial. Displacements have been 
scaled by a factor of 2 x 10 6 for ease of visualization. 

Figure [7] shows the u\ displacement component along a segment starting outside the core, to 
avoid the singularity that exists in the linearized elasticity solution. Note the close agreement 
between the analytic and numerical solutions for linearized elasticity. This is our first demon¬ 
stration of the viability of the proposed approach to model defects via dipole tensors, and exploit 
the native, variational structure of the theory of distributions. More such demonstrations follow 
in the coming sections. Also shown is the regularization of the displacement field obtained with 
the full, gradient theory of elasticity at finite strains. Note the significantly more gentle increase 
in \ui\ approaching the core. 

Figure [8] shows the solution, u\ through the core, highlighting the completely regularized 
gradient elastic solution at finite strain in comparison with the diverging, singular analytic 
solution. 

Figure [9] shows the trace of the first Piola-Kirchhoff stress tensor through the core of the 
point defect. Note the divergence of the classical, non-gradient, finite strain elasticity solution, 
in comparison with the regularization attained with even a very small gradient elastic length 
scale parameter, /. We also draw attention to the progressively lower variation in the stress field 
as the gradient elastic effect is strengthened. 

Figures an demonstrate the verstaility of our numerical approach to representing defect 
fields via dipole tensors in the variational, distributional setting. To the best of our knowledge, 
these figures also show the first numerical solutions to point defect fields with Toupin’s theory 
of gradient elasticity at finite strains. 
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Figure 6: The deformed configuration of a solid around an interstitial point defect. 
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Anal. soln. lin. elast. 
Num. soln. lin. elast. 
Num. soln. nonlin. elast. 


Figure 7: Point defect, u\ field along the segment from X = {0.55,0.5,0.5} to X = {1,0.5, 0.5}. 
The core is located at X = {0.5,0.5, 0.5}. The gradient elasticity result was obtained for finite 
strain with l — 0.05. 
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Figure 8: The u\ displacement field of a point defect along the segment from X = {0,0.5, 0.5} 
to X = {1,0.5, 0.5}. The analytical field diverges at the core. The gradient elasticity result was 
obtained for finite strain with l = 0.1. 
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Figure 9: The tr [P] stress field of a point defect along the segment from X = {0.4,0.4, 0.4} 
to X = {0.6, 0.6,0.6}. The classical, non-gradient, finite strain elasticity solution diverges, in 
comparison with the regularization obtained with gradient elasticity at finite strain. 
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5.2 Comparison of analytic and numerical solutions for the edge dis¬ 
location 


Consider an edge dislocation in the unit cube fio = (0, l) 3 , with half plane F' a subset of the 
X -2 — X 3 plane. The dislocation line and core are aligned with e;>, he at X -2 = 0.5 and have 
Burgers vector b = be\. We recall the analytic solution in the regime of linearized elasticity for 


such an edge dislocation in an isotropic, infinite medium (Hirth and Lothe 1982): 


uf 1 


_ _fr_r, , AW ¥ 2 

2t r [ Xi 2(1 - v)(Xl + X\y 


b 1 - 2 v 
27T 4(1 — v) 
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u% n = 0 


(45) 


The dipole tensor obtained by comparing our approach with Volterra’s solution (see Equation 
(23])) is 

2 /lb + Xb 0 0 

0 Xb 0 
0 0 Xb 


D = 


We apply the analytic displacement field ( |45| ) on the boundary surfaces T 0 as Dirichlet condi- 

(46) 


tions. The weak form (|30]) is re-written as: 

wi.JvudV - J wi,jDijdS = 0 


/ 


u = u 


on T n 


(47) 


Figure [lO] shows the distortion around such an edge dislocation, with displacements scaled by a 
factor of 20. 

Figure [TT] shows the trace of the u\ displacement component along a segment oriented with 
ei, but positioned slightly away from the core to avoid the singularity that exists in the Volterra 
solution (45). The Burgers vector has magnitude b = 0.008. The numerical solution with 


linearized elasticity was computed with linear, ( 7 °, basis functions and a knot span h = 0.016. 
We draw attention to the close match between this numerical solution and the analytic solution. 
Also shown is a trace of the u\ field computed with gradient elasticity at finite strain, using 
quadratic, C 1 , basis functions and the same knot span, h = 0.016. The gradient elastic length 
scale is l = 0.01. Note the regularization of the displacement field, which is the anticipated 
result of using a higher-order theory. 

Figure [12] shows the Pn component of the Piola-Kirchhoff stress tensor computed with 
gradient elasticity at finite strain as the length scale parameter is varied. Note the increasing 
degree of regularization of solutions as l increases from zero. For Z = 0, the classical, non¬ 
gradient, finite strain elasticity solution is obtained, and the stress diverges at the core. This 
singularity is reflected in the sharp increase in stress magnitude approaching the left end of the 
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Figure 10: The deformed configuration of an isotropic, linear elastic solid around an edge disloca- 








Figure 11: The u i displacement field for the edge dislocation along the line segment from X = 
{0, 0.51, 0.5} toX = {1,0.51,0.5}. The gradient elasticity result was obtained for finite strain with 
l = 0.1. 
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interval. These solutions were obtained with the Dirichlet boundary condition u = 0 on the 
boundary X 2 = 0, the Burgers vector magnitude b = 0.01, and for quadratic, C 1 basis functions 
with knot span h such that b/h = 11.5. 



Figure 12: Line plots of the P\\ stress component along a line segment from X = {0.53, 0.52, 0.5} 
to X — {1,0.52,0.5}, as the gradient length scale parameter, l is varied. 
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5.3 Comparison of analytic and numerical solutions for the screw dis¬ 
location 


Consider a screw dislocation in the unit cube Do = (0, l) 3 , with half plane T' a subset of the 
X 2 — X 3 plane. The dislocation line and its core are aligned with the e 3 direction, lie at X 2 = 0.5, 
and the Burgers vector is b = be 3. We recall the analytic solution in the regime of linearized 


elasticity for such a screw dislocation in an isotropic, infinite medium (Hirth and Lothe 1982): 


= 0 


uT = 0 


Uo 11 = —- tan 




(48) 


27r X\ 

The dipole tensor obtained by comparing our approach with Volterra’s solution (see Equation 


(23)) is 


D = 


0 0 fib 

0 0 0 

fib 0 0 


(49) 


tions. The weak form (30) is re-written as: 


We apply the analytic displacement field (|48j) on the boundary surfaces T 0 as Dirichlet condi- 

(50) 


/ wi^jcrijdV - / (w 1 j 3 D 13 + w 3 ^D 31 )dS = 0 
Jn 0 Jr f 


u = u 


on Tn 


(51) 


Figure [13] shows the distortion around such a screw dislocation, with displacements scaled 
by a factor of 20 . 

Figure [14] shows the trace of the u 3 displacement component along a segment oriented with 
ei, but positioned slightly away from the core to avoid the singularity that exists in the Volterra 
solution (48). The Burgers vector has magnitude b = 0.008. The numerical solution with 


linearized elasticity was computed with linear, C° basis functions and a knot span h = 0.016. 
We draw attention to the close match between this numerical solution and the analytic solution. 
Also shown is a trace of the u 3 field computed with gradient elasticity at finite strain, using 
quadratic, C 1 basis functions and the same knot span, h = 0.016. The gradient elastic length 
scale l = 0.01. Note the regularization of the displacement field, which is the anticipated result 
of using a higher-order theory. 

Figure [15] shows traces of the P 2 3 component of the Piola-Kirchhoff stress tensor computed 
with gradient elasticity at finite strain as the length scale parameter is varied. Note the increasing 
degree of regularization of solutions as l increases from zero. For l = 0, the classical, non¬ 
gradient, finite strain elasticity solution is obtained, and the stress diverges at the core. This 
singularity is reflected in the stress magnitude approaching the left end of the interval. These 
solutions were obtained with the Dirichlet boundary condition it = 0 on the boundary X 2 = 0, 
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Figure 13: The deformed configuration of an isotropic, linear elastic solid around a screw dislocation. 
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Figure 14: The displacement field for the screw dislocation along the line segment from X = 
{0, 0.51, 0.5} to X = {1,0.51,0.5}. The gradient elasticity result was obtained for finite strain with 
l = 0 . 01 . 
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the Burgers vector magnitude b = 0.01, and with quadratic, C 1 basis functions with knot span 
h such that b//h = 11.5. 



Figure 15: Line plots of the P 23 stress component along a line segment from X = {0,0.5, 0.5} to 
X — {1,0.5, 0.5}, as the gradient length scale parameter, l is varied. 


5.4 Energy of a single edge dislocation 


From classical, linearized elasticity we have the expression for the self energy of a single edge 
dislocation that is coaxial with a cylinder of radius R and incorporates a core cutoff, ro, to 
eliminate the singularity (Hirth and Lothe 1982): 


W„ 


fib 2 , R 

—7-7 l n — 

47t(1 — is) ro 


(52) 


Usually, this expression is presented as a logarithmic divergence with R. Instead, we have 
computed the self-energy of a single edge dislocation in the unit cube and carried out mesh 
refinement studies. As the knot span shrinks, quadrature points move closer to the core, and 
the computed self energy can be expected to diverge for fixed b as h 0 in the same manner 
that \n(R/ro) diverges for fixed R as ro 0. The boundary conditions for this computation, 
and for the remaining energy studies are U\ = 0 on X 2 = 0, with homogeneous traction and 
higher-order traction on the remaining boundaries. 

Figures [16] and [lT| respectively, show the strain energy and strain gradient energies of a 
single edge dislocation computed with gradient elasticity at finite strains. These computations 
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use b = 0.01. Note that for the gradient length scale l —» 0, i.e. for large values of 5/Z, the strain 
energy diverges logarithmically with b/h. This represents the regime wherein gradient effects are 
vanishingly present, and the classical, non-gradient response may be expected. Interestingly, the 
strain gradient energy also displays the same logarithmic behavior with b/h. For larger values 
of /, i.e. for b/l —>• 0 the strain gradients strongly regularize the problem and the logarithmic 
divergence with b/h is mollified for the strain energy as well as the strain gradient energy. We 
note that as l increases, both components of the elastic free energy decrease. This is due to the 
increasing stiffness of the displacement fields as l increases, which also is reflected in the above 
field solutions around defects. 


5.5 Studies of the elastic free energy of pairs of interacting, parallel 
dislocations 

In a cylindrical domain of radius R , the interaction energy per unit length between pairs of 
parallel dislocations, each with Burgers vector b and radial separation r, when calculated using 
classical, linearized elasticity, is 
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(55) 


(56) 


where a is a numerical parameter controlling the size of the core cutoff. Figure [18] compares the 
above interaction energy versus r for a pair of oppositely signed edge dislocations with computa¬ 
tions of the total energy of a pair of similarly interacting edge dislocations using our formulation 


for classical, linearized elasticity. Shown are the analytic result (53), the numerically computed 


result, and a shifted numerical result that accounts for the fact that the numerical domain is a 
unit cube with a square cylindrical core cutoff in comparison with the analytic solution, which 
uses infinitely long circular cylinders for the domain and the core. This shifting is achieved by 
requiring the solutions to match for large separations between the dislocations. Without this 
shifting, although the numerical solution, based on classical, linearized elasticity, has the correct 
trend, the values are systematically in error because of the difference in representation of the 
domain and core, and the combination of interaction and self energies in the numerical solution. 
We note that the energy of interacting dislocations decreases as the fields of the oppositely 
signed dislocations compensate with decreasing separation r. 

Figures [19}|2T] show corresponding results for like signed edge dislocations, oppositely signed 
screw dislocations and for like signed screw dislocations, respectively. In each case the trend 
shown by the numerical solution is correct, and the values improve upon the shifting as explained 
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Figure 16: The strain energy of a single dislocation, computed with gradient elasticity at finite 
strain plotted on (a) a linear scale for b/h , and (b) a logarithmic scale for b/h to emphasize the 
logarithmic divergence in the l —>> 0 regime (large b/l ). Note the convergence with mesh refinement, 
and the decrease in strain energy with an increase of gradient length scale, l (small b/l). 
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Figure 17: The strain gradient energy of a single dislocation, computed with gradient elasticity at 
finite strain plotted on (a) a linear scale for b/h , and (b) a logarithmic scale for b/h to emphasize the 
logarithmic divergence in the l —>> 0 regime (large b/l ). Note the convergence with mesh refinement, 
and the decrease in strain energy with an increase of gradient length scale, l (small b/l). 
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Figure 18: Energy of two edge dislocations with opposite signs; classical, linearized elasticity. 
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above to account for the domain and core shapes, and combination of interaction and self energies 
in the numerical solutions. In general, interacting oppositely signed dislocations show a decrease 
in total energy while like signed dislocations show an increase in energy. The results in Figures 
have been generated with b = 0.008 and a knot span h = 0.016 for linear, C° basis 
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functions. The analytic comparisons are with Equations (53 56) using a = 1. 
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Figure 19: Energy of two edge dislocations with like signs; classical, linearized elasticity. 
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Figure 20: Energy of two screw dislocations with opposite signs; classical, linearized elasticity. 
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Figure 21: Energy of two screw dislocations with like signs; classical, linearized elasticity. 
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Finally, Figures [22}|25] are equivalent computations for pairs of interacting, oppositely signed 
edge, like signed edge, oppositely signed screw and like signed screw dislocations, respectively, 
computed with gradient elasticity at finite strains, using b = 0.01 and b/h = 11.5. Here also, 
rather than attempt to define an interaction energy for nonlinear, finite strain elasticity, we 
simply present the total strain energy of the configuration versus the dislocation separation. No 
core cutoff is necessary in these computations because of the regularization of the singularity 
by gradient elasticity. Note the increased suppression of variation of the energy as l increases, 
which is a result of the regularization. 


x io' 3 




Figure 22: Energy of a pair of interacting edge dislocations with opposite signs; gradient elasticity 
at finite strain. 


5.6 Edge dislocation near a free surface 

Consider an edge dislocation in the unit cube flo = (0, l) 3 , with half plane T' a subset of the 
X 2 — Xs plane. The dislocation line and core are aligned with e 3 , lie at X\ = 0.95, X 2 = 0.5 
and have Burgers vector b = be 1 , where b = 0.01. Dirichlet boundary conditions u = 0 were 
applied at X 2 = 0, and the remaining surfaces are traction free. The schematic is shown in 
Figure [26j Our studies show convergence of the u\ displacement field with mesh refinement for 
different gradient length scales (Figure |27|), and of the Pn stress component (Figure [28]) as well. 
For both fields, u\ and Pn, an increased regularization is observed for larger /. For the Burgers 
vector used here, a knot span h > 0.1 fails to resolve the force dipole in the core. Note that Pn 
does not constitute the full traction component along ei, and therefore does not vanish at the 
traction free edge. 
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Figure 23: Energy of a pair of interacting edge dislocations with like signs; gradient elasticity at 
finite strain. 
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Figure 24: Energy of a pair of interacting screw dislocations with opposite signs; gradient elasticity 
at finite strain. 
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Figure 25: Energy of a pair of interacting screw dislocations with like signs; gradient elasticity at 
finite strain. 



es 

Figure 26: Schematic of an edge dislocation near a free surface, showing the half plane and force 
dipole. 
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Figure 27: Displacement component u\ plotted along the line segment between X — {0,0.5, 0.5} 
and X — {1,0.5, 0.5}. Note the convergence with h and the increasing regularization with /. 


5.7 The elastic field of a nonplanar dislocation loop 


The schematic of a non-planar dislocation loop in a finite domain is shown in Figure [29] with 
the corresponding half-planes of force dipoles. The Burgers vector of the dislocation loop is 
b = be i, where b = 0.001. Each linear segment has length 0.2, symmetrically located about 
the point X = {0.5, 0.5, 0.5}, with edge and screw segments discernible by the line direction 
relative to the Burgers vector. Dirichlet boundary conditions u = 0 are applied at X 2 = 0, and 
the remaining surfaces are traction free. 

Figure |30] shows contours of the P 12 stress component created by the non-planar dislocation 
loop. Note the higher stresses induced for smaller values of gradient length scale /, due to the 
loss of regularization. This example is motivated by a similar computation by |Roy and Acharya] 


(2005) using a continuum theory based on the dislocation density tensor. 


5.8 Low angle grain boundary represented by a distribution of edge 
dislocations 

Figure [31] is a schematic of a low angle grain boundary with a distribution of three edge dis¬ 
locations to represent the misorientation. The dislocations are parallel, along e 3 , with cores 
located at {Xi,X 2 } = {0.51,0.3}, {0.5, 0.7} and {0.49,0}. The Burgers vector is b = be 1 , with 
b = 0,0025. The gradient length scale is l = 0.0005. Dirichlet boundary conditions u = 0 are 
applied at X 2 =0, and the remaining surfaces are traction free. The grain misorientation is 
shown in Figure [32] and the distortion around the grain boundary in Figure |33| 
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Figure 28: Stress component P\\ plotted along the line segment between X — {0,0.5, 0.5} and 
X = {1,0.5, 0.5}. Note the convergence with h and the increasing regularization with /. 
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Figure 29: Schematic of a non-planar dislocation loop with (a) the Burgers vector, edge and screw 
segments and (b) corresponding half planes with force dipoles. The two vertical half planes are 
semi-infinite. 

6 Discussion and conclusions 

The crux of this work is that variationally based numerical methods are well-suited to repre¬ 
senting the elastic fields of defects. The key, enabling insight is that the force dipoles used 
classically to model defect cores give rise to singular distributions, which can only be manip¬ 
ulated in the variational setting; i.e., in weak form. Our numerical results match well with 
analytic expressions for defect fields in the limit of classical (non-gradient) linearized elasticity, 
when the dipole tensor is suitably parameterized. The extension to nonlinear elasticity poses 
no special difficulty, and when gradient elasticity is included, the expected regularization of the 
elastic fields is obtained. This cluster of results establishes the viability of representing defect 
fields in a fully numerical framework by modelling their cores as dipoles within distributional 
theory, and exploiting the variational framework. Beyond this, the interaction energy calcula¬ 
tions indicate that the approach has potential for modelling defect interactions. Finally, the 
free surface, dislocation loop and grain boundary examples suggest that computations on more 
complex defect configurations will also prove successful. 

It is important to note that by equating the core force dipole tensor with the effective force 
dipole tensor that emerges from linearized elasticity for point defects and dislocations, we have, 
in effect, matched those classical solutions. For the point defect, the comparison was made 
against the relaxation volume tensor, which is defined in the far-field. However, for dislocations, 
the dipole tensor was matched to the effective body force of Volterra’s solution. Nevertheless, 
these are not the best fits, since the linearized elasticity solutions are poor representations in 
the core. In ongoing work we are exploring the extraction of the dipole tensor by comparison 
with density functional theory calculations of defect configurations. Another front on which 
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Figure 30: Contours of P \2 induced by the non-planar dislocation loop. 


43 
























Figure 31: (a) A distribution of dislocations modelling a low angle grain boundary, (b) Half planes 
(red) corresponding to the three edge dislocations. 



Figure 32: The u\ field in neighboring grains. The deformation has been scaled by a factor of 20 
to make the grain misorientation apparent. Locations 1-4 are magnified in Figure [33| 
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Figure 33: The mesh lines demonstrate the distortion around the grain boundary. The deformation 
has been scaled by a factor of 10 for ease of visualization. The four sub-figures correspond to the 
locations in Figure [32] 
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we see potential for our approach is in computations involving large numbers of dislocations. 
The full elastic fields, including the influence of dislocation interactions will result from a single 
numerical solve without the need for multipole expansions. This could offset the computational 
expense of mesh resolution around defects that our method does require. Our method also works 
seamlessly to represent the interactions of point and line defects such as considered by ?. Even 
more attractively, these computations will be with gradient elasticity at finite strains. 

The evolution of large ensembles of dislocations is, of course, central to an understand¬ 
ing of mesoscopic material response. Our approach can be extended in this direction: The 
Peach-Koehler on a dislocation can be obtained via a variational calculation and emerges as a 
configurational force. The approach is similar to Steinmann (2002), with the additional feature 
of gradient elastic contributions. A formulation for dislocation dynamics would then be at hand 
with our formulation if the velocity of each dislocation were to be written via a first-order rate 
law driven by the Peach-Koehler force. The expense of multipole expansions could be avoided 
because a single numerical solve of the nonlinear elasticity problem would provide the driving 
forces on all dislocations, as already mentioned in the previous paragraph. 

There remain at least a few open questions in the context of our approach: (a) Can a 
variational calculation demonstrate a dipole tensor representing a single dislocation splitting 
into multiple dipoles, each representing partials? (b) The dipole tensor representation excludes 
any explicit representation of the separation of delta functions. For this reason, another open 
question is how we may model core spreading, such as in Zhang et ah (2015). One possible 
approach would be to switch to an explicit representation of Dirac delta monopoles with a finite 
separation, and carry out a variational calculation seeking to minimize the free energy with 
respect to the separation as a configurational variable, (c) As suggested by a reviewer, there 
remains the problem of whether a solid passing from a stress-free state to a stressed state upon 
nucleation of a dislocation represents a global minimum of total free energy, or whether the 
non-convexity of realistic free energy functions implies that the dislocated solid is in a local, 
not global, minimum. In this regard, we have carried out calculations based on earlier work 
(Rudraraju et al. 2014), which raise the possibility that non-convex free energy functions may 
allow such local, but not global minima of dislocated crystals. Another problem of interest is the 
interaction of defects with grain and phase boundaries. A continuum defect-based treatment 


was recently presented for this problem using dislocation and disclination densities (Acharya 


and Fressegeas, 2015). A straightforward extension of our treatment would be possible in which 


the energy stored in grain boundary configurations such as those in Section |5~8] would interact 
with neighboring defects, thus repelling or attracting them. A similar treatment could also be 
established for diffuse interface models of phase boundaries based on the approach in |Rudraraju| 


et al. (2015), where strain-derived order parameters establish the change in crystal structure 


between phases. The elastic energy interactions between these interfaces and defects modelled 
as shown here could be a viable representation of coupling of defects with phase boundaries. 
The sharp interface problem would require different approaches to include interface energies. 

We note that finite element computations of a full, field theory of the continuum dislocations 
have appeared in Roy and Acharya (2005), including comparisons of stress fields and a number 
of dislocation configurations. Their formulation consists of the governing equilibrium equations, 
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a div-curl system for elastic incompatibility arising from the definition of the dislocation density 
tensor, and a first-order wave propagation equation for the evolution of the dislocation den¬ 
sity. The incompatibility and wave propagation equations are solved with a least-squares finite 
element approach. Our formulation, based on the distributional dipole tensor and its rigor¬ 
ous variational treatment, may be considered somewhat simpler since a standard finite element 
method will do—at least for classical, non-gradient elasticity—with extra quadrature points. It 
is for the C 1 requirement of gradient elasticity that we invoke isogeometric analysis, although 
we have used it also for classical, non-gradient elasticity. 
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